Superfluid to Mott-insulator transition in Bose-Hubbard models 
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We study the superfluid- insulator transition in Bose-Hubbard models in one-, two-, and three- 
dimensional cubic lattices by means of a recently proposed variational wave function. In one di- 
mension, the variational results agree with the expected Berezinskii-Kosterlitz-Thouless scenario of 
the interaction-driven Mott transition. In two and three dimensions, we find evidences that, across 
the transition, most of the spectral weight is concentrated at high energies, suggestive of pre-formed 
Mott-Hubbard side-bands. This result is compatible with the experimental data by Stoferle et al. 
[Phys. Rev. Lett. 92, 130403 (2004)]. 
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Recent experiments on cold atoms trapped in optical 
lattices demonstrated that the Mott transition (MIT), 
originally introduced in electronic systems, [l[ can be ex- 
perimentally realized also in bosonic systems, f3| where 
the MIT is actually a superfluid-insulator transition. Re- 
cent experiments by Stoferle et al. [3] have shown that 
a considerable amount of spectral weight is concentrated 
at high energy even within the superfluid phase. More 
specifically, the data suggest that, especially in three- 
dimensions, a Mott-Hubbard gap of order U develops 
already on the superfluid side of the MIT, akin to what 
is predicted to occur in electronic systems. Q Although 
these evidences are not incompatible with the accepted 
theory of the critical behavior across the superfluid-to- 
insulator transition, 5] they clearly demand for a more 
detailed comprehension that must include also high- 
energy excitations. There have been already several the- 
oretical attempts, mainly based on suitable extensions of 
mean-field theory, to uncover the whole dynamical be- 
havior across the MIT. 0, 0] These calculations predict 
in the most general cases the existence of high energy 
modes even in the superfluid phase that might explain 
the experimental data. However, when the transition is 
approached at fixed integer filling by tuning for instance 
the interaction strength, these theories also predict that 
all modes soften at the transition. In this work, we intend 
to address this question by an alternative approach based 
on a variational wave function that has been recently 
proposed in the context of the electronic MIT. [lO| 
The accuracy of the wave function is checked by compar- 
ison with Green's Function Monte Carlo (GFMC) simu- 
lations, that allow us to obtain numerically exact results 
by a stochastic sampling of the ground-state wave func- 
tion. [ll| In contrast with the aforementioned mean-field 
theories, we find that, at fixed density, the MIT is accom- 
panied by a gradual transfer of spectral weight from the 
low-energy sound mode towards high energies, so that, 
when the Mott insulating phase is established, most of 



the spectral weight is already concentrated at high en- 
ergy. In addition, our analysis uncovers features of vari- 
ational wave functions able to describe a Mott transition 
that are novel and might be common to bosonic as well 
as fermionic systems. 

Bosons in optical lattices can be modeled by the Hub- 
bard Hamiltonian: 0, [H, [H, [H | 



U ^ , 
+ "2 2^ 



1), (1) 



where a] 



r) creates (annihilates) a particle at site 



with integer spin a = —S, ■ ■ . ,S, and = ^la'^itT- 
The Hubbard model H]) at integer filling has generally 
two different phases: one superfluid, for U < Uc, and the 
other insulating above Uc- In this work, we shall focus 
on the spinless case and attempt to describe the MIT 
by means of a variational wave function. In spite of the 
fact that the variational approach is a simple and well 
established technique, its application to the MIT turns 
out to be extremely difficult. For instance, the celebrated 
Gutzwiller wave function is not appropriate to describe 
the MIT, as it leads to an unrealistic insulator with no 
density fluctuations. 1^ lB| 



Recently, an extension of the Gutzwiller wave func- 
tion has been proposed [§|, that proved to be very ac- 
curate to describe an electronic MIT in ID. [^, [l^l Here 
we apply the same variational approach to the S = 
Bose-Hubbard model ([!]) with nearest-neighbor hopping 
t/2 in a one-dimensional chain (ID), a two-dimensional 
(2D) square lattice and a three-dimensional (3D) cubic 
lattice with L sites and periodic boundary conditions. 
We consider the following ansatz for the variational wave 
function 

I*) = exp 1^-^ X] Vijn^nj + gMB ^ 6 j l*o), (2) 
where j^o) is the non- interacting fully-condensed wave 
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FIG. 1: Variational results for the Jastrow potential Vq mul- 
tiplied by in ID and 2D and by \q\^ in 3D for increasing 
values oiU/t (from bottom to top). Upper panel: ID case for 
60 and 100 sites. Middle panel: 2D case for 20 x 20, 26 x 26, 
and 30 x 30 clusters (along the (1,0) direction). Lower panel: 
3D case for 8 x 8 x 8, 10 x 10 x 10, and 12 x 12 x 12 clusters 
(along the (1, 0, 0) direction). 
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FIG. 2: Density structure factor A'^, divided by \q\ calculated 
with variational Monte Carlo (left panels) and GFMC (right 
panels) in ID (upper panels) and 2D (lower panels) In ID, 
L = 60 (full symbols) and L = 150 (empty symbols) and 
Ujt = 1.6, 1.8, 2, 2.2, 2.4, 2.5 and 3. In 2D, L = 12 x 12 (full 
symbols) and 16 x 16 (empty symbols) for the variational 
calculation {U jt = 10, 10.2, 10.4, 10.6, and 10.8) and for the 
GFMC calculation {U jt = 8, 8.2, 8.4, 8.6, and 8.8). AU cases 
are shown from top to bottom for increasing values of I] jt. 



function, i.e. |$o) = (^I=o)^l^)' being h\ the creation 
operator at momentum k and N = L the number of 
particles. The components of the Jastrow potential, 
Vij = v{\Ri — Rj\), are independently optimized by a 
Variational Monte Carlo (VMC) minimization of the to- 
tal energy. 17| In the following, we will denote by pq and 



Indeed, since at least at weak-coupling, the expression 



Vq the Fourier transforms of the boson-density rii and of 
the Jastrow parameters Vij, respectively. Finally, qmb is 
a variational parameter related to the many-body opera- 
tor = hi Y{g{l - di+s) + di Y{s{l - hi+s), where hi = I 
[di = 1) if the site i is empty (doubly occupied) and oth- 
erwise, and 5 is the vector that connects nearest-neighbor 
sites. [l8| This term is kept just to improve the variational 
accuracy (mainly in 2D and 3D) but does not introduce 
important correlation effects, that are instead contained 
only in the long-range tail of the two-body Jastrow po- 
tential Vij . Both the Jastrow factor and the many-body 
operator ^ commute with the particle number, hence 
belongs to the Fock space with N = L bosons. 

Although our ultimate scope is to uncover some dy- 
namical properties across the MIT, we think it is worth 
discussing in some detail the role of the Jastrow factor 
in ([2]), which is the novel ingredient with respect to the 
conventional Gutzwiller wave function. In the superfluid 
phase, a long-range Jastrow potential is surely needed to 
restore the correct small-q behavior of the static density 
structure factor, i.e., Nq = - \q\. 



l+lVqN° 



(3) 



holds with 7 = 2, [19| and because the non-interacting 

~ i^olP-qPql^o) COnSt, it follows that Vq ~ 

Assuming that the expression ([3]) remains valid 
for \q\ — > even in the Mott insulating phase, one might 
be tempted to beheve that Vq ~ is necessary and 
sufficient to recover in any dimension the appropriate 
Nq insulating behavior, consequence of exponen- 

tially decaying correlation functions. However, one easily 
realizes that, were this conclusion correct, the variational 
wave function ^ could not describe any bosonic insula- 
tor in 3D, since Vq ^ is not sufficient to empty the 



condensate fraction. [20| Indeed, as shown below, the op- 
timized variational wave function has a more diverging 
Vq ~ l/I^P in the 3D Mott insulator, although Nq ^ q'^, 
implying that the formula ^ does not generally hold. 

In Fig. [1] we draw the optimized Jastrow potential Vq. 
For any dimension, the MIT is clearly signaled by the 
sudden change in the small-g behavior of Vq. On the 
one hand, the superfluid phase is always described by 
Vq ^ ct/\q\, with a increasing with U . On the other hand, 
the Mott insulator has a much more diverging Vq. In ID 
we recover the Vq ~ l/?^ behavior, like in the fermionic 
case. Q In 2D, the leading behavior of the Jastrow po- 
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FIG. 3: Upper panel: Density structure factor Nq calcu- 
lated by the variational Monte Carlo for 3D and U/t — 20. 
Lower panels: Nq for non-optimized wave functions with 
Vq ~ /Ssu/lgl^ for two values of /^sd. 



tential across the transition is less clearcut than in ID. 
Indeed, we cannot establish whether, on the insulating 
side, the leading behavior is given by Vq ~ (^20/0^ with 
/32D large but finite, or logarithmic corrections have to 
be considered, i.e., Vq ^ \n{l / \q\) / q'^ . Finally, in 3D a 
more diverging Vq l/l'Zp is stabilized in the insulating 
regime. Therefore, in all cases the Jastrow potential is 
able to empty the condensate. 2]J We note that, within 
this approach, the MIT shows up in the wavefunction in 
the form of a binding-unbinding transition of opposite- 
charged particles (empty and doubly occupied sites). Q 

In order to check the validity of our approach, we com- 
pare the VMC results of the density structure factor Nq 
with the numerically exact ones obtained by GFMC. At 
small q's we can generally write Nq ~ 7i|?| + 72'?^- In 
the superfluid phase, 71 ^ while, in the Mott insu- 
lator, 71 = and 72 7^ 0, see Fig. [2] In ID, we have 
evidence that 71 has a very sharp crossover from a finite 
value to zero across the MIT, suggestive of a true jump 
in the thermodynamic limit. Moreover, our numerical re- 
sults indicate that 72 diverges as the MIT is approached 
from the insulating side (this is particularly evident from 
the GFMC results). Within the variational approach, 
this behavior follows from Vq ~ Pm I q^ in the insulating 



phase with — > at the transition. In conclusion, 
the ID MIT can be located at Udt ^ 2.45 ± 0.05 in the 
VMC, whereas GFMC gives Uc/t ~ 2.1 ± 0.1 (in agree- 
ment with previous calculations of Ref. [Hjllii), showing 
that the wave function 12]) is not only qualitatively but 
also quantitatively correct. 

The density structure factor Nq displays quite distinct 



FIG. 4: Variational results for the average energy of density 
excitations E{q) of Eq. (gl in ID, 2D, and 3D. 



long-wavelength behaviors for weak and strong interac- 
tion also in 2D, see Fig. [2l The VMC structure fac- 
tor goes like Nq ^ "yi\q\ + 729^, for U/t < 10.3, while 
above we find Nq ~ 729^- The critical value of the 
on-site interaction is slightly different from the GFMC 
one, Uc/t ~ 8.5, which agrees with Ref. 2j|. In spite 
of slightly different values of Uc, the qualitative behavior 



across the MIT is similar both in VMC and in GFMC. 
Differently from ID, approaching the MIT in 2D, 71 — > 
while 72 7^ is smooth across the transition. 

Even more interesting is the 3D case. Here, the GFMC 
is severely limited by small sizes and, therefore, we just 
discuss the variational results. As we mentioned, the 
optimal Jastrow potential turns from Vq ~ ct/\q\ in the 
superfluid phase, into Vq ~ /^sd/I^P in the Mott insu- 
lator (see Fig. [1]). The sudden change of behavior al- 
lows to locate the transition around Uc/t ~ 18, which 
is close to the critical value of recent Monte Carlo sim- 
ulations in 3D. [2H1 Even though Vq l/I^P, the struc- 
ture factor in the Mott insulator has the correct behav- 
ior Nq ^ q^. In turns this implies that Eq. ([3]) does not 
hold, not even for \q\ 0, which is quite unexpected. 
In order to prove more firmly that Vq ^ /Ssd/I^P can 
indeed lead to Nq ~ g^, we have calculated the latter 
with a non-optimized wave function of the form ([2]) with 
Vq ~ /?3D/|9p and for different values of /Jan- As shown 
in Fig. [21 for small we have Nq ~ \q\^, implying that 
Eq. ([3]) is qualitatively correct. However, above a crit- 



ical P^jj, the behavior turns into Nq 



q , signaling a 



remarkable breakdown of Eq. ([3]). The optimal value of 
that we get variationally at the MIT is larger than 
P^jj, confirming our variational finding Nq ~ q^. We note 
that the change of behavior as a function of P^d is consis- 



4 



tent with the binding-unbinding phase transition recently 
uncovered in a classical 3D gas with l/jgl '^-potential. [26[ 
Let us now come back to our original motivation con- 
cerning the dynamical properties across the superfluid- 
insulator transition. Although our variational wave func- 
tion is meant only to describe the ground state, it can also 
provide important insights into the structure of the exci- 
tation spectrum. One can easily prove that the following 
expression holds in model H]): 

D ,n s I T\ dujuj S{q,uj) 
i=i ^ duj S{q,Lj) 

where (T) is the average value of the hopping, the sum 
is over the spatial directions, D = 1, 2, 3 is the space di- 
mensionality, and S{q,uj) the dynamical structure factor. 
E{q) is the first moment of S{q, uj) and can be regarded 
as the average energy of density excitations, which is 
therefore directly accessible through our variational cal- 
culation. In the superfluid phase E{q) oc q at small q, 
while E{q) develops a finite gap, i.e., E{0) ^ 0, in the 
Mott insulator, which is an upper bound of the actual 
Mott-Hubbard gap. In ID this gap seems to vanish at 
the MIT, in agreement with the Berezinskii-Kosterlitz- 
Thouless scenario, see Fig.lH On the contrary, both in 2D 
and 3D, we find that E{0) is finite and of order U every- 
where in the Mott insulating side, even right at the MIT, 
see Fig. |4l This implies that high-energy excitations ex- 
ist in the Mott insulator and carry most of the spec- 
tral weight. Also interesting is the behavior of the linear 
slope of E{q) within the superfluid phase as the MIT is 
approached. We recall that, assuming for the structure 
factor the small-g expression Ng =71 \q\ + q"^ + 0{q^), 
71 is finite at the MIT in ID, while it vanishes in 2D and 
3D. Moreover, as the MIT is approached, 72 diverges 
in ID but stays finite in 2D and 3D. Since the hopping 
energy is finite and continuous across the transition, it 
follows that in 2D and 3D the linear slope of E{q) should 
diverge at the MIT, although our numerical evidence is 
more clearcut in 3D than in 2D. Excluding the possibil- 
ity that the sound velocity diverges at the MIT, we must 
conclude that the spectral weight is gradually transferred 
from the sound mode to high-energy excitations that ex- 
ist already in the superfluid phase and are smooth across 
the MIT, suggestive of pre- formed Mott-Hubbard side 
bands. These results are actually consistent with the ex- 
perimental data of Ref . [sj . 

In conclusion we have demonstrated that a long-range 
Jastrow potential does allow for a faithful variational de- 
scription of a Mott transition in the bosonic Hubbard 
model. The average energy of the charge-density exci- 
tations, that is accessible by our calculation, suggests 
in two and three dimensions that pre-formed Hubbard 
side-bands coexist with sound modes in the superfluid 



phase near the Mott transition and carry most of the 
spectral weight in the insulator. This is an interesting 
and also surprising result, that bears a lot of similarities 
with the MIT in electronic systems, 0] but is not ac- 
counted for by most accepted theories of the superfluid 
to Mott insulator transition in bosonic systems. In anal- 
ogy with the bosonic example, we expect that a singular 
Jastrow potential Vq ~ is necessary to describe 

the 3D Mott transition in fermionic models, too, all the 
more reason when realistic Coulomb interaction is taken 
into account, in which case a Jastrow potential 1/q^ is 
necessary already in the metal to reproduce the proper 
long-wavelength behavior. 
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